!
!  DARKBAT     Dark frame subtract utility.
!  Copyright (C) 1997-x  Filip Hroch, Masaryk University, Brno, CZ
!
!  This file is part of Munipack.
!
!  Munipack is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, either version 3 of the License, or
!  (at your option) any later version.
!  
!  Munipack is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!  
!  You should have received a copy of the GNU General Public License
!  along with Munipack.  If not, see <http://www.gnu.org/licenses/>.

program darkbat

  use input
  
  implicit none

  ! No. of image dimensions
  integer, parameter :: DIM = 2
  
  ! Version:
  Character(len=*),parameter :: version = 'DARKBAT Version 2.0, (C)1997-2010 F.Hroch (hroch@physics.muni.cz)'

  ! The darkbat assumes that output image has in intensity
  ! within range minvalue < x < maxvalue.
  ! The minvalue always is 0 for measured intensities (?)
  ! The maxvalue is determined as (2**bitpix - 1).
  ! The specific characteristic of dark images may require some modification
  ! for real-type intensity of image.
  real, parameter :: minvalue = 0.0

  character(len=32768) :: sciname,darkname,outname,buf,key,value, & !mask,
       line,record
!  character(len=15) :: form
!  integer :: icit = 1
!  integer :: i,j, n, nj, mj, nhead
  integer :: i,j, n, nhead, io
  real :: maxvalue
  real, dimension(:,:), allocatable :: DARK, SCI, OUT
  character(len=80), dimension(:), allocatable :: head

! FITS variables
  integer :: blocksize,naxis,naxes(DIM),pcount,gcount,bitpix,bpix,bitpix0
  logical :: simple,extend,anyf!,exist

  ! input of top line with init values
!  read(*,*,end=90,iostat=j) darkname,mask,icit,bitpix
!  if( j /= 0 ) stop 'Bad format of input data.'

  bitpix0 = 0
  i = 0

  do
     read(*,'(a)',iostat=io,end=20) record
     if( io /= 0 ) stop 'Failed to read of an input.'
  
     call input_rcomment(record,line)

     select case(input_classify(line))

     case(input_key)

        call input_keys(line,key,value)

!        if( key == 'Output mask' ) then
!           read(value,*,iostat=i) mask
!
           ! prepare output mask
!           mj = index(mask,"?")
!           if ( mj > 0 )then
!              do nj = 0, len_trim(mask)
!                 if( mask(mj+nj:mj+nj) /= '?' ) exit
!              enddo
!              write(form,"(A,I1.1,A,I1.1,A)") "(I",nj,".",nj,")"
!           endif
!           outname = mask
!
!        endif

        if( key == 'Darkframe' ) then
           read(value,*,iostat=io) darkname
           
           ! open darkframe
           call ftopen(20,darkname,0,blocksize,i)
           call ftghpr(20,DIM,simple,bpix,naxis,naxes,pcount,gcount,extend,i)
           if( i /= 0 ) goto 666
           allocate(DARK(naxes(1),naxes(2)),stat=j)
           if( j /= 0 ) stop 'Not enought memory for dark.'
           call ftg2de(20,1,minvalue,naxes(1),naxes(1),naxes(2),DARK,anyf,i)
           call ftclos(20,i)  
           if( i /= 0 ) goto 666

           ! allocate buffers for scientific exposure
           allocate(SCI(naxes(1),naxes(2)),OUT(naxes(1),naxes(2)),stat=j)
           if( j /= 0 ) stop 'Not enought memory for scientific image.'

        endif
             
        if( key == 'Output BITPIX' ) then
           read(value,*,iostat=io) bitpix0
        endif

        if( key == '' .or. value == '' .or. io /= 0 ) then
           write(*,*) 'This record: ',trim(key),trim(value),' parse failed.'
        endif


     case (input_data)
        
        read(line,*,iostat=io) sciname,outname
        if( io /= 0 ) then
           write(*,*) 'Failed to read of filename:',trim(line)
           stop 'Bad format of the input filename.'
        end if


!  do
     ! read name of scientific exposure to subtract
!     read(*,*,end=90,iostat=j) sciname
!     if( j /= 0 ) stop 'Bad format of input filename.'
        write(*,'(5a)',advance="no") '(',trim(sciname),') - (',trim(darkname),') = '

        ! generate output filename:
        ! by mask
!        if( mj > 0 ) then
!           write(outname(mj:mj+nj-1),form,iostat=i) icit
!           if( i /= 0 .or. index(outname,'*') > 0 ) &
!                stop 'The output filename mask doesn not conforms to image number order.'
!           icit = icit + 1
!        endif

        ! temporary
!        if( mask == '.' ) then
!           outname = sciname
!        endif

        ! work in current dir
!        if( mask == '$' ) then
!           j = index(sciname,'/',.true.) + 1
!           outname = sciname(j:)
!        end if

        if( .not. allocated(DARK) ) &
             stop 'Dark frame not loaded yet.'

        ! scientific exposure initialisation
        call ftopen(25,sciname,0,blocksize,i)
        call ftghpr(25,DIM,simple,bpix,naxis,naxes,pcount,gcount,extend,i)
        if( i /= 0 ) goto 666
        if( naxes(1) /= size(DARK,1) .or. naxes(2) /= size(DARK,2) )then
           write(*,*) "The image and darkframe size does not match. Skipping."
           call ftclos(25,i)
           goto 665
        endif
        
        call ftghps(25,nhead,n,i)
        allocate(head(nhead),stat=j)
        if( j /= 0 ) stop 'Not enought memory for head.'
        do j = 1, nhead
           call ftgrec(25,j,head(j),i)
        enddo

        call ftg2de(25,1,minvalue,naxes(1),naxes(1),naxes(2),SCI,anyf,i)
        call ftclos(25,i)
        if( i /= 0 ) goto 666

        ! dark subtraction
        OUT = SCI - DARK

        if( bitpix0 == 0 ) then
           bitpix = bpix
        else
           bitpix = bitpix0
        end if

        ! range checking
        if( bitpix > 0 )then
           maxvalue = 2**bitpix - 1
           where( OUT < minvalue )
              OUT = minvalue
           endwhere
           where( OUT > maxvalue )
              OUT = maxvalue
           endwhere
        end if

        ! check if output file exist, does not overwrite
!        Inquire(file=outname,exist=exist,iostat=i)
!        if( exist )then
!           ! remove file to overwrite
!           if( mask == '.' .or. &
!                ( mask(1:4) == 'dout' .and. mask(10:13) == 'fits' ) ) then
!              open(26,file=outname,iostat=j)
!              close(26,status='delete',iostat=j)
!              if( j /= 0 ) then
!                 write(*,*)"Error occured during remove of output file. Skipping."
!                 goto 665
!              endif
!           else
!              write(*,*) "The output file exist (",trim(outname),"). Skipping."
!              goto 665
!           end if
!        endif

        
        ! write out corrected image
        call ftinit(26,outname,blocksize,i)
        do j = 1, nhead
           call ftprec(26,head(j),i)
        enddo
!        if( bitpix /= bpix ) then
        call ftukyj(26,'BITPIX',bitpix,'number of bits per data pixel',i)
!           if( bitpix < 0 ) then
!              call ftdkey(26,'BSCALE',i)
!              call ftdkey(26,'BZERO',i)
!           else
!        endif
        if( bitpix > 0 ) then
           call ftukyj(26,'BSCALE',1,'',i)
           call ftukyj(26,'BZERO',2**(bitpix-1),'',i)
        endif
        call ftphis(26,'Darkframe subtracted: '//trim(darkname),i)
        call FTPCOM(26,version,i)
        call ftp2de(26,1,naxes(1),naxes(1),naxes(2),OUT,i)
        call ftclos(26,i)
        if( i /= 0 ) goto 666
        
        write(*,*) '(',trim(outname),')'
        
665     continue
        if( allocated(head) ) deallocate(head)

     end select

  end do

!     enddo
!90   continue
20 continue

!  deallocate(DARK,SCI,OUT)

666 continue

  if( allocated(DARK) ) deallocate(DARK)
  if( allocated(SCI) ) deallocate(SCI)
  if( allocated(OUT) ) deallocate(OUT)

  if( i /= 0 ) then
     write(*,*) 'The following error(s) has been occured:'
     do
        Call ftgmsg(buf)
        if( buf == ' ' ) exit
        write(*,*) trim(buf)
     enddo
  endif

end program darkbat
